测绘学报 | 邓开元：多光谱光学遥感影像水体提取模型

Original 测绘学报智绘科服 2022-07-16

收录于合集 #《测绘学报》 340个

本文内容来源于《测绘学报》2021年第10期（审图号GS（2021）6119号）

多光谱光学遥感影像水体提取模型

邓开元, 任超

桂林理工大学测绘地理信息学院, 广西桂林 541004

基金项目：国家自然科学基金(42064003)

摘要：地表水监测是遥感的重要基础应用。光学遥感水体提取的原理是基于不同的地物具有不同的光谱反射特性，某些地物（冰雪、阴影、云）因为与水体具有相似的反射特性导致提取分类失败。本文针对传统水体指数在水体提取时出现错分、漏分的问题，提出一种归一化多波段水体指数NDMBWI。分别用3个试验对本文指数的稳定性进行测试。试验1的区域为西藏林芝地区，数据源为同一时相的Landsat 8、Sentinel 2卫星影像，试验结果验证了本文指数抑制冰雪的能力，本文指数的Kappa系数为0.86、总体精度为0.93、错分误差为0.03、漏分误差为0.12、制图精度为0.97、生产者精度为0.88，均优于已有的指数。试验2的数据源为高分一号，以香港迪士尼为试验区域，在存在少量云的环境下进行水体提取，证明了本文指数能抑制云及其云下阴影。试验3提取了多个地区的水体，验证了本文水体指数的稳定性。本文使用多源光学遥感影像验证了NDMBWI的可行性，不需额外借助辅助数据，即可排除雪、云、阴影的影响，能够更加有效地自动化提取水体，可推广到海岸带资源研究、冰川变化、内陆湖泊变化等领域。

关键词：遥感光学遥感影像水体指数水资源水体提取

引文格式：邓开元, 任超. 多光谱光学遥感影像水体提取模型[J]. 测绘学报，2021，50(10)：1370-1379. DOI: 10.11947/j.AGCS.2021.20200482DENG Kaiyuan, REN Chao. Water extraction model of multispectral optical remote sensing image[J]. Acta Geodaetica et Cartographica Sinica, 2021, 50(10): 1370-1379. DOI: 10.11947/j.AGCS.2021.20200482

阅读全文：http://xb.sinomaps.com/article/2021/1001-1595/2021-10-1370.htm

引言

地表水变化对生态、社会和经济的影响是当前研究的热点^[1^]。光学遥感水资源研究主要可以分为：气候和水资源^[2^]、不同类型水体的提取^[3^-5^]、洪水监测^[6^]、河流形态^[7^]、水体长时间序列^[8^]、海岸带变化^[9^]。目前，Landsat影像是用于变化分析，识别和检测地表水的最广泛使用的数据源^[10^]。在基于Landsat遥感影像水体提取的研究中，多数是针对大区域的河流、湖泊等水体的提取^[11^-14^]。近年来，一些学者采用Sentinel 2影像对大区域的地表水进行提取研究，取得了良好的效果^[15^-16^]。目视解译是最初的光学遥感影像提取水体的方法^[17^]，其缺点是效率低、成本高，没有利用不同地物具有不同反射率的特性。遥感影像的水体信息提取是基于电磁波谱上的反映特性^[18^]，文献[19]通过使用Landsat 4影像的6个波段提出缨帽变换湿度指数(tasseled cap wetness, TCW)。受到归一化差异植被指数^[20^]的启发，文献[21]基于Landsat影像的绿波段和近红外波段提出了归一化差异水体指数(normalized difference water index, NDWI)，NDWI是简单的水体指数, 该指数能较好地抑制植被信息，但受薄云和山体阴影影响较大。文献[22]使用短波近红外波段代替近红外波段，改良了NDWI，提出了改进的归一化差异水体指数(modified normalized difference water index, MNDWI)。MNDWI指数能突出水体信息，但也容易把水体和阴影混淆。基于分类二叉树的思想，文献[23]提出了针对有阴影的场景自动水体提取指数AWEIsh和针对没有阴影的场景自动水体提取指数AWEInsh，二者都是多波段水体指数。AWEI指数在提出时，声明该指数不适用于高反射率的地物，比如冰雪、白色建筑物等。文献[24]在全球大范围多个地区进行测试了自身的稳定性和可行性, 提出了水体指数(water index 2015, WI2015)。针对复杂环境，文献[25]成功地构建了一个多波段水体指(multi-band water index, MBWI)，MBWI一共使用了5个波段并且在多个复杂环境的地区测试了MBWI的有效性。MBWI能较好地抑制非水体信息的能力，但其缺点之一是不适用提取细小水体。面向对象提取水体信息的方法^[26^], 在抑制阴影信息的同时，能够准确提取细小水体。但是面向对象的方法的缺点之一是：不具有普遍性，难以在大范围下实现高精度的水体信息提取。

在地形复杂的环境下，如高山积雪较多的区域，对水体提取产生干扰的因素有：阴影、冰雪等^[27^-28^]。如何减少遥感影像的噪声干扰仍是水体信息提取的难点。本文分析多种地物的反射率，综合考虑多波段间的关系，通过建立六元不等式组，设计一个顾及多个地物噪声影响的光学遥感影像归一化多波段水体指数(normalized difference multi-band water index, NDMBWI)。

1 多光谱影像水体指数模型构建

1.1 已有的水体指数

选取7种典型的、常用的水体指数进行研究，7种常用的水体指数见表 1。这7种水体指数以绿波段和近红外波段、短波近红外波段为主要影响成分。TCW、MNDWI、AWEI、WI2015主要考虑绿波段和短波近红外波段的反射率，因为水体的绿波段反射率大于短波近红外波段的反射率，而其他非水体地物的绿波段小于短波近红外波段的反射率。正常情况下，这些指数都可以区分水体与非水体。但是，TCW、MNDWI、AWEIsh、WI2015会把可见光波段反射率大于短波近红外波段的非水体地物分类为水。文献[25]指出，雪就是可见光波段的反射率远远大于短波近红外波段的地物。

表 1 已有的水体提取指数Tab. 1 Commonly used water extraction index

指数	来源	等式
TCW	文献[19]	0.031 5Blue+0.202 1Green+0.310 2Red+0.159 4Nir-0.680 6Swir1-0.610 9Swir2
NDWI	文献[21]
MNDWI	文献[22]
AWEInsh	文献[23]	4(Green-Swir1)-0.25(Nir+2.75Swir2)
AWEIsh	文献[23]	Blue+2.5Green-1.5(Nir+Swir1)-0.25Swir2
WI2015	文献[24]	1.720 4+171Green+3Red-70Nir-45Swir1-71Swir2
MBWI	文献[25]	2Green-Red-Nir-Swir1-Swir2

表选项

1.2 模型的机理

本文的技术路线如图 1所示。在不同的时段、不同的地区，从多幅Landsat 8、Sentinel 2、高分一号(GF-1)影像里选择不同地物不同波段的真实地表反射率光谱值进行拟合，具体采样点的信息见表 2。内陆水体、海洋水体、阴影、雪、植被、旱地、建筑物、云的平均反射率如图 2所示。图 2曲线是基于表 2中的采样点求出来的平均值。水体的反射主要在蓝绿光波段，其他波段的吸收率很强，尤其是在近红外波段、短波近红外波段有强烈的吸收带。研究已有的水体指数可以从各个光谱间的关系作为切入点，本文以不同波段与绿波段的比值，阐述已有水体指数模型的机理及其存在的缺陷，并建立本文水体指数NDMBWI。表 3为各个光谱与绿波段的关系，表 3中的数值可看成图 2曲线量化后的结果。图 1中，可以把TCW、NDWI、MNDWI、AWEIsh、AWEInsh、WI2015近似看成双波段水体指数。水体指数的模型机理在于不同波段的对不同地物具有不一样的敏感性和有效性。雪、云的反射特性与水体相似，表 3中，部分地物(内陆水体、海洋水体、雪、云)的Q₄、Q₅、Q₆的值都小于1，从机理上说明：主要依赖绿波段、近红外波段、短波近红外波段的双波段水体指数在有雪环境下不能有效地提取水体。

图 1 本文技术路线Fig. 1 Technical route of this paper

图选项

表 2 数据采样点涉及的数据Tab. 2 Data involved in data sampling points

卫星类型	地区	时段	采样点的数量	涉及的地物
Landsat 8	林芝	2020-02-13	50	内陆水体、雪、旱地
Sentinel 2	林芝	2020-02-13	50	内陆水体、雪、旱地
Sentinel 2	北京	2019-12-11	50	内陆水体、建筑物
Sentinel 2	巴彦淖尔	2019-09-27	50	内陆水体、旱地、植被
Sentinel 2	和田	2019-09-22	50	内陆水体、旱地、植被、云
Sentinel 2	青岛	2019-12-28	50	浑浊海洋水体(河口)、建筑物
Sentinel 2	齐齐哈尔	2019-09-04	50	浑浊内陆水体、建筑物
Sentinel 2	三亚	2019-11-24	50	海洋水体、植被、建筑物
Sentinel 2	梧州	2019-09-20	50	内陆水体、植被、建筑物
Sentinel 2	香格里拉	2019-03-29	50	内陆水体、雪、阴影、植被
GF-1	香港	2018-10-02	50	海洋水体、云、阴影

表选项

图 2 拟合多个地物的反射率曲线Fig. 2 Reflectance curves of multiple features

图选项

表 3 各个光谱关系Tab. 3 Relation of each spectrum

地物	蓝波段与绿波段的比值Q₁	绿波段自身的比值Q₂	红波段与绿波段的比值Q₃	近红外波段与绿波段的Q₄比值	短波近红外1波段与绿波段的比值Q₅	短波近红外2波段与绿波段的比值Q₆
内陆水体	0.74	1	0.94	0.44	0.18	0.49
海洋水体	1.33	1	0.77	0.64	0.6	0.49
阴影	0.87	1	1.02	1.35	1.22	0.71
雪	1	1	1	0.92	0.1	0.11
植被	0.45	1	0.64	6.83	3.54	1.67
旱地	0.65	1	1.3	1.64	2.1	1.86
建筑物	0.81	1	1.08	1.2	1.32	1.21
云	0.95	1	1.03	1.12	0.93	0.75

表选项

表 3中的数值是来源于Landsat 8、Sentinel 2、GF-1在各个地物、各个波段上的平均值，然后各自相除得到的比重。本文试验数据源为Landsat 8、Sentinel 2、GF-1，具体信息见表 4。

表 4 数据源信息Tab. 4 Data source information

波段	Landsat 8中心波长	Sentinel 2中心波长	GF-1中心波长
蓝波段	482.5	490	502
绿波段	562.5	560	576
红波段	655	665	680
近红外波段	865	842	810
短波近红外1波段	1610	1610	无
短波近红外2波段	2200	2190	无

表选项

表 4中，Landsat 8、Sentinel 2、GF-1三者的中心波长不完全一致，这也导致了已有的水体指数模型在用于不同的数据源时，表现出阈值跳动的问题。经过研究，已有的水体指数的模型可用式(1)表示

(1)

式中，Y表示水体指数的值；x₁、x₂、x₃、x₄、x₅、x₆是6个波段前的系数；ρ₁至ρ₆是蓝波段、绿波段、红波段、近红外波段、短波近红外1波段、短波近红外2波段的反射率。

式(1)整体除以绿波段，可变形为

(2)

根据水体的水体指数的值理论上要大于α，非水体的水体指数的值理论上要小于α。式(2)中，ρ₁、ρ₂、ρ₃、ρ₄、ρ₅、ρ₆的值用比值系数代替, 得到以下不等式组

(3)

式中，a、b、c、d、e、f、g、h分别代表求得的比值系数。

在本文把α设置为0，目的是希望本文指数能够以稳定的阈值提取水体信息，到达自动化提取水体的目的。结合式(3)与表 3的参数，组成一个六元不等式组

(4)

式中，x₁、x₂、x₃、x₄、x₅、x₆分别是蓝波段、绿波段、红波段、近红外波段、短波近红外1波段、短波近红外2波段的系数。

式(4)的系数源于Landsat 8、Sentinel 2、GF-1在各个地物、各个波段上的平均值，若是直接赋予式(4)每项等于定值，会因为推导不严谨，导致水体指数模型出现错漏现象，本文通过最小二乘法的方式进行消元、拟合，求出不同定值下的解。

由表 3中可以发现，已有水体指数Q₁、Q₂、Q₃、Q₄、Q₅、Q₆与1的大小关系并不能完全地实现水体的识别。例如，NDWI由绿波段与近红外波段组成，即是从Q₄的值可以解析NDWI对水体信息识别的能力，大于1的数值的地物包括阴影、植被、旱地、建筑物、云，小于1的值包括内陆水体、海洋水体、雪，此时导致NDWI无法以固定的阈值区分雪、水体。MNDWI、AWEIsh、AWEInsh、WI2015近似看成由绿波段与短波近红外波段组成的模型，可从Q₅的值可以解析此类型模型对水体信息识别的能力，大于1的数值的地物包括阴影、植被、旱地、建筑物，小于1的值包括内陆水体、海洋水体、雪、云。此时，以上指数只考虑Q₄或Q₅(双波段)，并没有综合考虑Q₁、Q₂、Q₃、Q₄、Q₅、Q₆的内在联系，所以本文建立六元不等式组，综合考虑多波段间的关系，设计一个顾及多个地物噪声影响的水体指数模型。

式(4)没有精确解，其原因是：式(4)是超定方程组的类型，即方程个数大于未知量个数。式(4)的解集是一个区域，解多元不等式(组)就是求它的各个不等式的解集的交集。为了使本文水体指数适用于国产遥感影像(没有短波近红外1波段、短波近红外2波段)，需要提前把x₅、x₆设置为0。本文通过逐步消元、拟合，求出其中的一种解：x₁、x₂、x₃、x₄、x₅、x₆分别是-1、3、2、-5、0、0，即初步构建的水体指数为

(5)

式中，ρ₁、ρ₂、ρ₃、ρ₄代表蓝波段、绿波段、红波段、近红外波段的反射率。

式(5)的系数是可行域中的简化解，目的是为了更好推广NDMBWI以及方便记忆。另外，可以通过使用式(5)的参数，反解式(4)进行验证。因为双波段水体指数难以抑制雪、云的影响，本文指数没有将绿波段作为唯一的正参数，额外引入了红波段作为正参数。为了使本文指数归一化，最终拟定的水体指数模型为

(6)

式中，ρ₁、ρ₂、ρ₃、ρ₄代表蓝波段、绿波段、红波段、近红外波段的反射率。

本文推算出来的系数，实则为式(6)的分子。在提取水体时，起决定性作用的是分子而不是分母。这样做的好处是：预处理时，在大气校正环节，某些算法为了减少存储空间会把结果扩大10 000倍。此时，本文使用了归一化的形式，则不会受以上的影响，无论大气校正有无扩大10 000倍，NDMBWI的范围始终在-1至1的范围。使用本文水体指数和表 1中列举的已有的水体指数进行仿真模拟，在不同的地区、不同的数据源选择水体、阴影、雪、植被、旱地、建筑物、云等7种地物作为对象，每种地物选择100个样本点，不同水体指数的稳定性测试结果如图 3所示。图 3表明，0可以作为NDMBWI的阈值，能够有效地区分水体与多种非水体；以上的水体指数均可区分水体与植被、旱地、建筑物；除了NDMBWI、MBWI，其他的水体指数均不能区分水体与冰雪；MBWI会对部分水体漏分；NDWI的阈值提高到0.2左右可以分区水体与冰雪。

图 3 不同水体指数的稳定性测试结果Fig. 3 Stability test results of different water indexes

图选项

2 水体提取实例分析

2.1 不同卫星影像在有雪环境下的水体提取

试验1选择的地区位于西藏林芝地区。数据源是Landsat 8影像和Sentinel 2影像，影像拍摄时间均为2020年2月13日。图 3中得知：8种指数以阈值为0均能将水体与部分非水体(植被、旱地、建筑物)区分。所以在本次试验中，NDMBWI、MBWI、WI2015、AWEInsh、AWEIsh、MNDWI、TCW统一以0作为阈值进行水体提取，而NDWI则是使用0.2作为阈值。提取雅鲁藏布江水体的结果如图 4所示。图 4(a)为Landsat 8影像的原始图。图 4(b)为Sentinel 2影像的原始图。对比图 4(c)、(d)可以发现，NDMBWI应用于Landsat 8影像、Sentinel 2影像时，均可以抑制了雪的影响，完成了提取水体信息。

图 4 雅鲁藏布江水体提取结果Fig. 4 Water extraction results of Yarlung Zangbo river

图选项

在试验1中，选择100个水体、非水体样本组成混淆矩阵，计算出的Kappa系数、总体精度、错分误差、漏分误差、用户精度、生产者精度见表 5，使用Sentinel 2影像的效果均优于Landsat 8。为了进一步研究水指数的阈值，研究NDMBWI、MBWI、WI2015应用于Sentinel 2影像的数值分布，三者的直方图如图 5所示。在直方图中，最优阈值一般选定在“双峰”之间的低谷区间。图 5(a)验证了在试验1中，0可以作为NDMBWI的最优阈值。图 5(c)并无明显的低谷，说明WI2015本身难以抑制冰雪的干扰。

表 5 分类精度混淆矩阵分析Tab. 5 Classification accuracy confusion matrix analysis

影像	指数	Kappa	总体精度	错分误差	漏分误差	用户精度	生产者精度
Landsat 8	NDMBWI	0.86	0.93	0.03	0.12	0.97	0.88
MBWI	0.62	0.81	0.16	0.24	0.84	0.76
WI2015	0.52	0.76	0.36	0.04	0.68	0.96
AWEIsh	0.58	0.79	0.25	0.14	0.75	0.86
AWEInsh	0.46	0.73	0.33	0.12	0.67	0.88
MNDWI	0.60	0.80	0.28	0.04	0.72	0.96
NDWI	0.66	0.83	0.24	0.04	0.76	0.96
TCW	0.42	0.71	0.30	0.29	0.70	0.71
Sentinel 2	NDMBWI	0.88	0.94	0.03	0.10	0.97	0.90
MBWI	0.76	0.88	0.05	0.20	0.95	0.80
WI2015	0.52	0.76	0.32	0.04	0.68	0.96
AWEIsh	0.54	0.77	0.30	0.06	0.70	0.94
AWEInsh	0.46	0.73	0.34	0.08	0.66	0.92
MNDWI	0.60	0.80	0.28	0.04	0.72	0.96
NDWI	0.68	0.84	0.24	0.02	0.76	0.98
TCW	-0.09	0.46	0.59	0.78	0.41	0.22

表选项

图 5 NDMBWI、MBWI、WI2015的直方图Fig. 5 Histogram of NDMBWI, MBWI and WI2015

图选项

2.2 有云环境下的高分一号影像水体提取

云一直是光学遥感影像难以处理的问题之一，本文提出的NDMBWI适用于存在少量云的情景。云一般都有阴影，传统的水体指数无法处理阴影。由仿真模拟的结果可知：NDMBWI能够以0为阈值区分水体与云、阴影。GF-1是中国自主研制的第一颗民用高分辨率光学卫星，为我国水环境的遥感监测提供更高空间分辨率的光学遥感数据源。但是，GF-1只有可见光波段与近红外波段，使用了短波近红外波段的MBWI、WI2015、AWEI、MNDWI、TCW等指数不能应用在GF-1。试验2以中国香港的迪士尼公园为试验区域，影像的拍摄日期为2018年10月2日。使用NDMBWI提取水体的对比结果如图 6所示。

图 6 香港迪士尼地区水体提取结果Fig. 6 Water extraction results in Hong Kong Disneyland

图选项

图 7作为NDMBWI在试验2的直方图，验证了0可以作为NDMBWI的最优阈值。云是光学影像无法处理的痛点，去云步骤一般分为两步：①识别云像素；②把云像素设置为背景值(无效值)。本文在设计指数模型时把云考虑到非水体地物的因子，所以通常情况下NDMBWI在会把云识别为非水体地物，云通过NDMBWI得到的值为负值。即NDMBWI无须做去云步骤。但是，值得说明的是NDMBWI适用于云量为20%以下的光学遥感影像。也就是说NDMBWI是无法适用于云量过多的场景。

图 7 NDMBWI的直方图Fig. 7 histogram of NDMBWI

图选项

2.3 多地区水体提取

选取表 2中的齐齐哈尔、巴彦淖尔、梧州、三亚作为试验3，从多个不同的地区，复杂的地形环境下提取不同类型的水体，验证本文水体指数模型的稳定性，图 8中的NDMBWI均以0作为阈值。

图 8 多地区水体提取结果Fig. 8 Water extraction results in multiple regions

图选项

3 结论和讨论

如何减少噪声干扰对遥感影像的不利影响是水体提取的难点，本文设计一种归一化多波段水体指数NDMBWI，避免积雪、阴影等的影响。试验1使用相同时相、同一地区的Landsat 8影像、Sentinel 2影像，以西藏林芝地区为试验区域，试验结果说明NDMBWI能抑制雪的影响。试验2使用GF-1影像，在影像含有少量云的情况下完成水体提取，验证NDMBWI可以抑制云、阴影的影响。阈值，对于水体指数是至关重要的。从水体指数的构建开始，要注意阈值的选取问题。试验3选择多个不同地区进行试验，以说明本文水体指数模型的稳健性。今后，构建新的光学水体指数时，应该考虑多个波段间的敏感性和有效性。本文的创新点在于：充分利用可见光波段、近红外波段的信息，不需额外借助辅助数据，即可排除雪、云、阴影的影响，完成水体信息的提取。

从理论上，任何的水体指数都无法保证提取所有的水体，本文提出水体指数虽然在某些情况下能得到较高的精度，并在冰雪、阴影、建筑、云等区分具有一定的优势，然而因为自然水体的复杂性，并不是任意一种水体指数能区分所有水体。也就是说，NDMBWI的最优阈值不能一直保持在0，例如在钱塘江悬浮泥沙较多的出海口地区提取水体时，NDMBWI的阈值应该适当降低。目前主流的水体指数的研究大多数是固定阈值。本文的下一步工作是：基于原有的水体指数模型，构建动态阈值的模型。例如，以NDVI作为NDMBWI的阈值，使用NDMBWI与NDVI模型进行双通道并行处理，在NDMBWI的矩阵中，大于NDVI的像素点可以归为水体，反之为非水体。本文方法稳健，可以作为基于深度学习的水体检测分类的基础，解决难以区分雪、水的问题。

作者简介
第一作者简介：邓开元(1994-), 男, 硕士生, 研究方向为遥感图像算法。E-mail: ytkz11@163.com
通信作者：任超, E-mail: renchao@glut.edu.cn

初审：张艳玲

复审：宋启凡
终审：金君

往期推荐

资讯

○ 中国科学院空天信息创新研究院城市环境遥感团队岗位聘用职工招聘
○ 中国工程院院士郭仁忠：城市如何数字化？既要有好“路”还要有好“车”
○ 中国科学家首获地球观测组织卓越个人奖
○ 关于诚邀参加中国测绘科学研究院2021年测绘科技成果推介会的函
○ 北京市测绘设计研究院2022年应届毕业生公开招聘
○ PPT | 李春来：基于高光谱红外成像的气体探测方法、现状与展望
○ 中国测绘 | 筹建北京冬奥赛道的测绘⼈

会议

○ 第一届中国生态系统遥感学术研讨会（三号通知）
○ 《测绘学报》2021年编委会会议在青岛召开 ○ 会议通知 | 2021中国地理信息科学理论与方法学术年会通知（第二号）
○ 关于第八届“测绘科学前沿技术论坛”再次延期召开的通知

《测绘学报》

○《测绘学报》2021年第10期目录

○ 测绘学报 | 宗敬文：测高重力异常中央区奇异积分数值求积公式

○ 测绘学报 | 左溪冰 : 高光谱影像小样本分类的图卷积网络方法